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A propagation-dispersion equation is derived for the first passage distribution function of a particle 
' moving on a substrate with time delays. The equation is obtained as the continuous limit of the 

first visit equation, an exact microscopic finite difference equation describing the motion of a particle 
on a lattice whose sites operate as time-delayers. The propagation-dispersion equation should be 
contrasted with the advection-diffusion equation (or the classical Fokker-Planck equation) as it 
describes a dispersion process in time (instead of diffusion in space) with a drift expressed by a 
propagation speed with non-zero bounded values. The temporal dispersion coefficient is shown to 
exhibit a form analogous to Taylor's dispersivity. Physical systems where the propagation-dispersion 
equation applies are discussed. 
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j3 ' I. INTRODUCTION 
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Often to describe the microscopic mechanism of a diffusion process, one considers a test particle executing a random 
walk on some substrate, and one writes a mean-field equation in terms of the probabilities that the particle performs 
elementary displacements in given or arbitrary directions. The question one then asks is: where will the particle be 
I . after some given time (in the long-time limit)? The answer is given by the distribution function F(r, t), the probability 
that, given the particle was initially at r = at t = 0, it will be at position r at time t (for t — > oo, that is for t large 
compared to the duration of an elementary displacement). F(r,t) is obtained as the solution to the Fokker-Planck 
equation for diffusion, and one finds that, in the the long-time limit, F(r,t) is Gaussianly distributed in space [[l]. 

When there is an interactive process between the particle and the substrate such that the particle undergoes directed 
motion subjected to time delays, it is interesting to view the motion in terms of first passages, and one expects that 
the long-time dynamics will be different from that described by the usual diffusion equation. We obtain indeed a 
. new equation for the long-time behavior of the first visit distribution function f(r, t) of a particle whose dynamics is 
£S| ' governed by a distribution of time delays. The main results in this paper are (i) the propagation- dispersion equation 
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where c is the propagation speed, and (ii) the expression for the time dispersion coefficient in terms of the covariance 
of the reciprocal velocity fluctuations 

where £ is a characteristic correlation length. We first give a heuristic derivation of the propagation-dispersion equation 
using multi-scale analysis which is then further substantiated by a more mathematically rigorous development. Using 
the latter method, we show that the above results generalize to inhomogeneous systems with c — ► c(r) and 7 — * 7(7*). 

A characteristic example where the equation applies is particle dispersion in a granular medium as studied ex- 
perimentally by Ippolito et al. pi, as we shall discuss below. On the other hand, there are prototypical abstract 
systems, for particle-substrate interactive dynamics, such as the automaton known as Langton's ant ||] and other 
simple automata (i[ [|, where particle-substrate interactions produce time delays in the dynamics. These automaton 
systems offer the advantage of explicit microscopic dynamics which can be solved exactly. For instance Grosfils, Boon, 
Cohen, and Bunimovich |5(] developed a one-dimensional automaton for which they provided a mathematical analysis 
also applicable to the two-dimensional triangular lattice. 
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The one-dimensional case is particularly simple to describe. The automaton universe is the one-dimensional lattice 
where, at each time step, a particle moves from site to site, in the direction given by an indicator. One may think of 
the indicator as a 'spin' (f or J.) defining the state of the site: when the particle arrives at a site with spin up, it moves 
to the next neighboring site in the direction of its incoming velocity vector, whereas its velocity is reversed if the spin is 
down. But the particle modifies the state of the visited site (T^=H) so that on its next visit, the particle is deflected 
in the direction opposite to the scattering direction of its former visit. With this specific microscopic dynamics, 
back-scattering produces time delay. Grosfils et al. derived the equations describing the microscopic dynamics of the 
particle on the one-dimensional lattice (and also on the triangular lattice) under the general condition that the spins 
at the initial time can be arbitrarily distributed. They proved that the particle will always go into a propagation 
phase, regardless of the initial distribution of spins; they also showed that the basic mechanism for propagation is a 
blocking process ||. 

It was shown by Boon || that the propagation equations in |5) are particular cases of a general equation describing 
the first visit of the particle to a new site in the propagation phase. This equation, obtained in the context of specific 
lattice dynamics, is a first passage equation which has more general applicability, and will be our starting point in 
the present work. 



II. FIRST VISIT EQUATION 



Consider a particle propagating in a D-dimensional channel; we project its motion onto a one-dimensional lattice 
- the propagation line - whose sites are labeled by integers I = 0, 1, 2, ... . The distance, in lattice units, between 
neighboring sites on the propagation line is denoted by Sr. We set the clock at to = when the particle is at site I = 
where it enters the propagation channel. Its trajectory will intercept successively sites I = 1, 2, 3, ... for the first time 
at times t\, t 2 , t 3 ... respectively. The Vs are integer multiples of the automaton time step St. While sites 0, 1, 2, 3, ... 
are equally spaced, the time differences between first visits, U+i — U, are (in general) not equally distributed. 

A given (arbitrary) spin configuration defines a set of first passage times ti, t2, ■■■ ■ Let i r be the random variable 
which corresponds to the number of steps required for the particle to reach position r = l r Sr for the first time. We 
define f{r,t) as the probability of finding the particle at position r for the first time at time t = t r = i r St; f(r,t) 
obeys the first visit equation, i.e. the finite difference equation H 

n 

f(r,t) = J2Pjf(r-Sr,t-r j ). (1) 

3=0 

Here pj is the probability that the particle propagates from r — Sr to r in Tj/St time steps, i.e. Tj is the time 
between two successive first visits on the propagation line. The sum is over all possible time delays, weighted by the 
probability pj with Sj=o Pj = For specific lattice dynamics, this is a condition for the existence of a blocking 
mechanism (described in j5|) responsible for propagation, and an explicit expression can then be given for Tj (see 
Appendix). Here we present a general derivation for which this specification is not necessary. 
An alternative equivalent formulation of Eq. (|I]) reads 

oo 

f(r,t)=jr i p j f(r-5r,t-j5t), (2) 

3=0 

which we shall use in Section [vj. The difference between the two formulations is that the structure of the distribution 
of the delays is contained either in the Tj's (Eq. ([!])) or in the pj's (Eq.(^|)). 

Equation (|l|) - or equivalently @ - expresses the probability that the particle be for the first time at position r at 
time t in terms of the probability that it was visiting site r — Sr at earlier time t — Tj, where j = 0, 1, 2, n. Given a 
site r on the propagation line, the particle will infallibly reach that site in the course of its displacements; the question 
then is: when will the particle be at position r when r is large? Since the particle executes tj St displacements to cover 
the distance r, the answer will be given in terms of the time distribution of the probability f(r, t) for large fixed value 
of r, i.e. for l r ^S> 1. This corresponds to taking the continuous limit of equation ([!]) in a physically consistent way 



that we establish in Section III and in a mathematically rigorous way as shown in Section VI. 

For one particular realization, the successive time delays are set by a given spatial configuration of the time delayers, 
and the time taken by the particle to perform a displacement from r — Sr to r depends on that configuration. For an 
ensemble of realizations, the distribution function of the time delays defines the average displacement time 

n n 

( T ) = v i T J = Pj 5t = (fj) St i (3) 

3=0 ]=0 
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and the variance 

(r 2 )-(r) 2 = |:£^ft-Eww] 2 |(ft) a 
[j=0 j=0 J 

= ({^)-{tf){5t)\ (4) 

where fij = tj / St is the number of time steps during the time delay Tj . The general condition on the pj distribution is 
that its moments be finite (for specific lattice dynamics, such as described in the Appendix, they are finite by virtue 
of the blocking mechanism discussed in B). Higher order moments (r a ) are defined similarly. 



III. MULTI-SCALE ANALYSIS 



Here we discuss how the continuous limit should be taken given that the system exhibits two time scales which 
correspond to (i) a propagation process characterized by the average time necessary to complete a finite number of 
displacements r/Sr 

E[U] = (M>rJ£ , (5) 

and (ii) the dispersion around this average value characterized by the variance 

Var[i r ] = «,i 2 ) - (tf){5tf^-. (6) 

For finite r, these are finite quantities. Correspondingly we define the following quantitites that will be used in the 
continuous limit 

c Or 

and 



7 = «M 2 > - (M) 2 ) ^ ■ (8) 



c (7^ 0) will be identified as the propagation speed (see Section IV) and 7 (> 0) will be identified as the dispersion 
coefficient (see Section [v|) . 

We want to compute the continuous limit of the first visit equation ([!]), and obtain a partial differential equation 
for f(r,t). The procedure must be performed in two successive steps according to the scale over which one wants to 
probe the process when r is large. This is analogous to multi-scale expansion in the derivation of the Fokker-Plank 
equation @ or of the Navier-Stokes equation B. Consider that, to measure first passages, we use a detector with 
tunable resolution. The first visit time t r goes like r (see (^)); so in order that the measure be performed with the same 
accuracy at any position r, we need a resolution such that t r /r has always the same order of magnitude. We can then 
measure first visit times at any position r with a resolution which is appropriate to evaluate the average first passage 
time, i.e. to measure the propagation speed c. However, when performing measurements over many realizations, the 
successive arrival times at position r (the fluctuations around the average time) will be poorly resolved (see insets 
in Fig.l). In order to measure the dispersion around ~E[t r ] with sufficient accuracy, the detector must be adjusted so 
that dispersion measurements at various positions r can be performed with the same resolution. Therefore we impose 
that the width AT = y/Vai[t r ] of the dispersion curve be measured with a resolution ST such that the curve always 
contains about the same number of points, i.e. AT /ST has the same order of magnitude for all positions. Since 
AT grows with the distance like y/r (see @), ST must also go ~ y/r in order to obtain an acurate measure of the 
dispersion 7 (see (||)). 

We may summarize by saying that the continuous limit must be taken in a sequential order, first to obtain the 
average first passage time (which gives a measure of the propagation velocity c), and second, to measure the first 
passage dispersion in the moving reference frame (i.e. around E[t r }). The first step will yield an Euler type equation 
and the second step the propagation-dispersion equation. Mathematically, these two steps materialize in the two 
successive orders, 0{t r ) and C(e 2 ), of the development in Section |l^, where e is a smallness parameter defined as 
St /t b s , the ratio of the microscopic time to the macroscopic observation time. The continuous limit corresponds 
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to e <C 1. It follows from the above discussion that there are two length scales: (i) the first one corresponds to 
propagation 

1 Sr 

t obs c = t obs — — = e~ <5r , (9) 

for which the measurement unit is e~ 1 Sr; (ii) the second length scale corresponds to dispersion 

= _ w ^ = ((mVMT 1 ^, (io) 

which is thus measured in units e~ 2 5r. It is therefore natural to introduce two space variables, r\ = er and r2 — e 2 r, 
and correspondingly two time variables, t\ = et and t% = e 2 t. Note that r\ and ri are merely two expressions of 
the space variable r with different scalings, and similarly for t\ and ti. Consequently the space and time derivative 
operators must be rescaled as 

d r - ed ri + e 2 d r2 ; d t - eft, + e 2 3 f2 (11) 

Accordingly f(r,t) is expanded as / = J*- -* + e/^ + C(e 2 ), where J*- -* is the distribution function in the absence of 
dispersionf/^ ' plays the same role as the local equilibrium distribution function in the derivation of the Navicr-Stokes 
equation gj). 

While in taking the continuous limit, there is a mathematical separation of scales (i.e. one probes either the 
propagation mode or the dispersion mode), in practice, there are systems with sufficiently large 7 (see Section VII 
and Appendix) where propagation and dispersion can be measured simultaneously, provided one uses a detector with 
high resolution. 



IV. PROPAGATION-DISPERSION EQUATION 

In order to compute the continuous limit of Eq. ([!]) , it is convenient to rewrite the equation as 

n 

f(r + 8r,t) = J2p j f(r,t-T j ) 1 (12) 



1=0 



to which we now apply the multi-scale expansion (in Section ^ we give an alternative derivation which is more 
mathematically rigorous). 



To first order, we obtain 



(0) 



(13) 



or, with (fl), 



drJ {Q) + -d t jW =0, (14) 



which is an Euler equation. 
To second order, we have 



0(e 2 ): d n f® + 9 ri /« + \5rd 2 r J^ = 

-M (9t2/ (0) +9ti/ (i)) + iip^/W, (15) 
or 2 or 

where we will inject the first order result (p"3]). This substitution can be performed either in the last term on the 
l.h.s. or in the last term on the r.h.s. with two different results of which only one can be correct. In fact there is 
no ambiguity as to which is the legal procedure: Eq. ( |l2|) is local in space (not in time) which indicates that in the 
continuous limit we must obtain apartial differential equation with a spatial derivative not exceeding first order. The 
correct substitution of (H) into dj|) 

gives 



drjW + = -\ (d t J (0) + d t J (1) ) + \{j^- 5 i) C 



/(°) , (16) 
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where the quantity in parantheses in the second term of the r.h.s. is equal to ((^t 2 ) — (n) 2 )(8t) 2 / 8r = 7. 
Summing Eqs.(|l4|) and (|l^) after multiplication by e and by e 2 respectively, we obtain 



ed,,(/ (0) + e / (1) ) + ^ /C°) = -I eftl (/(")+ e/ d)) - ^ dt2 f^) + l e ^ tJ i9). (17) 

Rearranging terms (incorporating terms of negligible higher order) and going back to the original variables yields 

d r f(r,t) + -d t f(r,t) = \ 1 d 2 f(r 1 t) , (18) 
c 2 

where terms of order e > 2 are omitted. With the initial condition that at the origin, say at r — 0, f(0, t) — S(t), the 
solution to Eq.([l8|) is 

f(r,t) = /l (7r) -^ exp L . (19) 

V 2n \ 2jr J 

These results are confirmed by the more mathematically rigorous derivation given in Section Note that in the 
case that all pj's are zero except one (uj = fi, e.g. for Langton's ant), then 7 = 0, and (|l8|) reduces to the Euler 
equation (^3) with c = 5r/(fj,St) (see 0). 

It is clear from (|l8|), that c is a propagation speed, and 7 is a transport coefficient expressing dispersion in time 
(instead of space like in the classical Fokker-Planck equation for diffusion). Equation (18) is the propagation-dispersion 
equation governing the first-passage distribution function of a propagating particle subject to time delays. Propagation 
is guaranteed because c is non-zero ({//) is finite); c has a finite maximum value for n = (i.e. j — and pa — 1), in 
which case 7 = (see (|J)), and there is propagation without dispersion. 

In Fig.l we show that the above analytical results are in perfect agreement with the numerical solution of the 
general equation, Eq.(|]). We used delay times equally distributed (Fig. la) and exponentially distributed (Fig. lb); 
explicit values of the quantities c and 7 are given in the figure caption. 



V. DISPERSION AND CORRELATIONS 

(i) The dispersion coefficient. It follows from Eqs.@ and (Q), that 7 is given by 

7 . «> - <'"> 2 (20) 
r 

which, for large r, is reminiscent of the classical expression for the diffusion coefficient: D = lim t ^ 0O (r 2 (t))/2t. 
Comparison of the two expressions shows interchange of space and time, a feature which is illustrated in Figs. 2 and 3. 
In Fig. 2, we show three typical runs on the 1-D spin-lattice by plotting the first passage time (minus the mean time 
of arrival (t r ) = r/c) as a function of distance r. One observes that [t— (t r )] as a function of r exhibits fluctuations of 
the same nature as those obtained when plotting the position of a random walker as a function of time. In Fig. 3 we 
present simulation data illustrating Eq.(|20|). From measurements performed in the 1-D lattice, the variance (t 2 r ) — (t r ) 2 
is seen to be a linear function of distance (see inset) with a slope equal to 7 (main frame) in the same way as the 
diffusion coefficient is obtained as the slope of the mean-square displacement versus time in the long-time limit of the 
classical random walk. 

(ii) The correlation function. From Eqs.(^) and ([?]), we have | = ^r-. Here t r can be expressed in terms of the local 
propagation velocity v(r), a fluctuating quantity with average value c. In fact it is the reciprocal local velocity which 
is physically relevant: it is the time taken by the particle to propagate from position r to r + Sr (divided by 8r). Then 
indeed 

(*r) = (f dr'J-jr) = f dr'(J-) = ~, (21) 







v(r') Jo v(r') c 



which is consistent with the definition of the propagation speed. 
Expressing the time t r in terms of the local velocity v(r), we have 

<a " <I *' ^FT [ *" ^> " I " L *" ( ^F) ^> ■ (22) 
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FIG. 1: Probability distribution f(r = 3 x 10 4 , t) based on general equation (|l). Here fAj = 1 + 2j, so c = [(1 + 2(j})] 1 . (a) 
time delays equally distributed for j = 0, 1, • • • , 9, with pj — 0.1; c = 0.1 and 7 = 33; half-width = y/2'yr ~ 1.41 x 10 3 . (b) 
time delays exponentially distributed: pj = C exp—[3j, with j = 0, 1, 9, j3 = 0.25, and C = Ej-o-?'] -1 = 1/45; c = 0.128 
and 7 = 52.7; half-width = <J2yr ~ 1.78 x 10 3 . The numerical simulation data and the analytical expression (Eq.(^9|); solid 
line, not visible) coincide perfectly. Insets: see text. In this and subsequent figures, time unit is the automaton time step. 



and 

lf 2 \ - (i„\ 2 = r At 1 r dr"t( 

^v(r') v J \v(r") v 
In terms of the reciprocal velocity fluctuations (5u _1 (r) = w _1 (r) — (v^ 1 ) = v~ x (r) — c -1 , (||j|) reads 



(t 2 r ) {U) 2 = r dr> ^ dr"{(-±- {-)) f - 7-777 - <-) )) ■ (23) 



(t 2 r )-{t r ) 2 = dr' <h"(5v- 1 (r f )5v- 1 (r f ')). (24) 



The dynamics of the propagating particle implies that the correlation function on the r.h.s. of ( p4j ) is ^-correlated, 
i.e. (Sv -1 (r 1 ) Sv^ 1 (r")) = <f>Q <5(y — V) with O = {(Sv _1 ) 2 ) = (-7) — -7, and where £ is the elementary correlation 
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FIG. 2: Displacement time (t — (t r ) = t — r/c) as a function of distance r for three different realizations in the 1-D lattice (see 
Appendix: q = .5, and j = 0, 1). Vertical axis: time in automaton time steps; horizontal axis: distance in lattice unit lengths. 



length. So it follows from (gg) and (|§ that 

7 = £((0-c- 2 ), (25) 

that is 7 is the covariance of the reciprocal velocity fluctuations multiplied by the correlation length (here equal to 
one lattice unit length). This result is analogous to Taylor's formula of hydrodynamic dispersivity which is expressed 
as the product of the covariance of the velocity fluctuations with a characteristic correlation time B. We call 7 
the temporal dispersion coefficient. Evidently, there is no dispersion (7 = 0) in the absence of velocity fluctuations 
(00 = 0). 

From these results, the explicit computation of the dispersion coefficient is straightforward. For instance for the 
one-dimensional spin-lattice (see Appendix; Fig. 4a), the local velocity is either 1, with probability q, or 1/3, with 
probability (1 — q), and the reciprocal mean velocity is c _1 = q + 3(1 — q); so 

7 - q + 3 2 (l-q)-c- 2 = 4(7(1-9), (26) 

which is exactly the value of 7 evaluated for the 1-D spin lattice in the Appendix and in . A similar computation 
for the triangular lattice (see Fig. 4b) yields the value 7 = 72g(l — q) in accordance with the result obtained in (|]. 

(iii) The time current. Writing Eq.(|l8|) as 

{Or + -d t )f(r,t) + d tJ (r,t) = 0, (27) 
c 

with j(r, f) = — 9r d t /(r, t), shows formal analogy with the classical mass conservation equation 

(d t + c-V)p(r,t) + V-j(r,t) = 0, (28) 
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FIG. 3: Dispersion coefficient, see Eq.(j2Cj), in 1— D lattice (see Appendix: q — 0.5, j = 0, f , po = q, pi = f — q, and fij = l + 2j). 
Inset shows Iff) — (t r ) 2 as a function of r, for one single realization, and main frame shows 7 measured from the slope in inset; 
the horizontal dashed line is the theoretical value: 7 = ({(J 2 ) — (a*) 2 ) — 9+ (1 — 9)3 2 — (3— 2q) 2 — 4g(l — q) = 1. r is in lattice 
unit lengths and 7 is in automaton time steps squared per lattice unit length. 



with space and time variables interchanged. So in (p7[), j(r,t) can be interpreted as a "current in time ". 

(iv) The control parameter. In classical advection-diffusion phenomena, the control parameter is the Peclet number 
P = UL/2D, where U denotes the mean advection speed, L, the characteristic macroscopic length, and D, the diffusion 

) . The analogue for propagation-dispersion follows by casting Eq. ([l^) in non-dimensional form 

2T 

dr f(r, t) + d x /(r, t) = B- 1 d 2 t f(r, t) ; B = — . (29) 

7c 

Here r and t are the dimensionless space and time variables: r = r(cT) -1 and t = tT~ 1 , where T is a characteristic 
macroscopic time. B is the control parameter for propagation-dispersion: it is a measure of the relative importance 
of propagation with respect to dispersion. Indeed, B = — — ^—-L — Lp/Lp, i.e. the ratio of the characteristic 

"yc "y cl 

dispersion length Lrj to the characteristic propagation length Lp. At high values of B, i.e. Lp> 3> Lp, the distribution 
function is very narrow, and transport over large distances (r > cT) is dominated by propagation. 

(v) The power spectrum. The propagation-dispersion equation (Eq), subject to the initial condition f(r — 0, t) = 5(t), 



coefficient (see e.g. 10 
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describes an initial value problem with initial value fixed in space; so we can solve the equation by spatial Laplace 
transformation. Using k as the conjugate space variable, we obtain 

nf(K,w) - f(r = 0,u) = -—f( K ,w) - liu 2 f( K ,u), (30) 

c z 

where f(u>) denotes the time Fourier transform. With the initial condition f(r = 0, t) = S(t), i.e. f(r = 0,u) = j-, 
© yields 

] M = ±-{k + ^ + |. 2 )- 1 . (31) 

This result shows that the system dynamics is described by one single mode: k = — — — Z id 2 . The corresponding 
spectrum S(k,u>) = 2Ref(n = ik,uj) is given by 

S{KUJ) = ~L (k + u/c)"+ i( 7 ^)2 ' ( 32 ) 

with J" °° S(k,u>) dk = 1. The spectrum ( |32| ) (which can also be obtained by double Fourier transformation of (|lS|)) 
exhibits a single Lorentzian line typical of a diffusive phenomenon, but there are two essential differences with the 
spectrum obtained from the classical advection-diffusion equation: (i) the spectrum is a Lorentzian in k (rather than 
in w) with half-width at half-height Afc = ^joj 2 , indicating that dispersion is diffusive in time (instead of space); (ii) 
the Lorentzian is shifted by a quantity proportional to the reciprocal of the propagation speed. 

VI. EXACT SOLUTION OF THE FIRST- VISIT EQUATION 

In order to emphasize the discrete nature of the problem, we introduce the notation 

g(l,i) = (St) f(l5r,i6t), (33) 
where here, and below, latin arguments (Z, • • •) always indicate integers. The first visit equation, Eq.(|^), is then 

oo 

^M)=Ep?" x) $(J-M-i). ( 34 ) 

3=0 

Here we consider the general case where the transition probabilities p^p depend on the lattice position I as indicated 
by the superscript (for simplicity we omit the tilde notation). It is convienient to introduce the generating function 
for the distribution which is defined to be 

oo 

ht(x) = jV0(M). (35) 

i=0 

and from which the distribution is obtained via 

1 cf 

g{l,i) = lim -—hi(x). (36) 
x^o i\ ax 

Temporal moments of the distribution can be calculated as 

so that Mi(z) = hi(e z ) is the generating function for moments of the first passage time at lattice position I, The 
boundary condition that the particle starts at lattice site I = at time i = implies that <?(0, i) — dio, or equivalcntly 
h (x) = 1. 

Substituting Eq.@ into Eq.@ gives 

OO CO 



i = j = 

OO 

= hi-tix^p^aP, (38) 
j=o 
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so that the general solution is 



/ / oc 

hi( X ) = h Q (x) n £ 

k=l \j=0 



(fe-1) , 

Pj x J 



(39) 



Using this solution and Eq.(|37|), exact moments may be easily calculated; with the normalization J^°L Pj = 1, we 
obtain 



where 



(A*>J = 



1 / oo 

tEE 



I 

fc=l j=0 

1 1 

fe=l 



0-1) . 



E fc-1 .J 



E^ 1} i 



(40) 
(41) 



(42) 



(43) 



are just averages over the elementary process. In order to develop the limiting form of the distribution for large I, 
we assume that and of are of order 1 for all I as is certainly true if the elementary probabilities are independent 
of lattice position. We then introduce a new stochastic variable which measures deviations away from the expected 
waiting time as 



Wl 



u - mo 



flof 



The moment generating function Ni(z) for this new variable is related to that for the original variable by 

1((jl)iZ 



Ni(z) = exp 

or, using ( |39| ) with the boundary condition ho(x) = 1 

l{p)iz 



(44) 



(45) 



'l k=l \j=0 

By double expansion of the second term on the r.h.s. of ([46|), we obtain 

1 



lnN t (z) 



-z 2 + 



1 

71 



(46) 



(47) 



So in the limit of large I — r/Sr, the generating function for the moments of wi is just Ni(z) = exp (z 2 /2) which 
is recognized as the generating function for a Gaussian distributed variable with unit variance. We conclude that in 
this limit the distribution for the original (temporal) variable becomes 



5(U)z»i 



: exp 



i (i- mo 



/lof 



(48) 



from which the distribution for the physical quantities (for r 3> 8r) reads 



(r(r) r) 2 exp 



2T(r); 



(49) 
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with 



1 



r / Sr 



C(r) 
r(r) 



1 . . St 1 



r (Wr) = of 



a (**) 



2 



fe=l 

r / 5' 



Sr 



(50) 
(51) 



fe=i 



where the definition of (i)fe-i and (t 2 )k-i follows straightforwardly from (|42| ) and ([l3[). 

This result, ([l9]), is just a realization of the central limit theorem. Notice that the development is completely 
independent of any assumption on the magnitude of Sr and St. To make contact with the dynamical formulation in 
Section^ for the case of position-independent probabilities, we note that the present result is derived with the initial 
condition g(0, j) — Sjo, which means that it is the Green's function for the difference equation. In the continuous limit, 
the corresponding boundary condition is /(0, t) = S(t) so that the dynamical equation for the distribution in this limit 
should be linear and should have JIs| ) as its Green's function which implies the generalized propagation- dispersion 
equation 



with 



1 _ (t)r _ \ -> (r) • 

c(r) Sr ^— ' 3 Sr ' 

.7 = 1 



and 



7(r) = l(t 2 ) r (t)t 



E^HE^ 



(52) 



(53) 



Sf_ 

Sr 



(54) 



Notice the difference between c(r) and y(r) in (|52j) and C{r) and Y{r) in (49): I.e. symbols denote local quantities 
whereas capitals indicate space-averaged quantities. When the waiting time probabilities are space- independent, 
c(r) — > c and 7(7*) — > 7, and the generalized equation (52) reduces to Eq.(18). An example of microscopic dynamics 



whose continuous limit is described by the generalized equation is discussed in the Appendix. 



VII. COMMENTS 

There is an algebraic similarity in the structure of the propagation-dispersion equation ( |l8| ) and of the classical 
advection-diffusion equation jjj which can be formally transformed into each other by interchanging space and time 
variables. It should be clear that the two equations describe different, but complementary aspects of the dynamics 
of a moving particle. Solving the propagation-dispersion equation answers the question of the time of arrival and 
of the time distribution around the average arrival time in a propagation process. It is also legitimate to ask the 
complementary question "where should we expect to find the particle after some given time ?" which should be long 
compared to the elementary time step, but short with respect to the average time of arrival. We will then observe 
spatial dispersion around some average position which can be evaluated from the solution of the advection-diffusion 
equation. This observation stresses the complementarity of the two equations. 

Because the propagation-dispersion equation describes the space-time behavior of the first passage distribution 
function f(r,t), i.e. the probability that a particle be for the first time at some position, it describes transport where 
a first passage mechanism plays an important role. So the equation should be applicable to the class of front-type 
propagation phenomena where any location ahead of the front will necessarily be visited, the question being: when 
will a given point be reached? 

A most interesting case is the "Diffusion of a single particle in a 3D random packing of spheres" to quote the title of 
an article by Ippolito et al. Q where the authors describe an experimental study of the motion of a particle through 
an idealized granular medium. They measure particulate transport and 'dispersivity' which corresponds precisely 
to the quantity 7 computed in the present paper. In particular the experimental data presented by Ippolito et al. 
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show that the mean square transit time of the particle through the medium is a linear function of the mean transit 
time (Figs. 10 and 11 in B) itself a linear function of the percolating distance (Fig. 2 in this observation is an 
experimental illustration of the features described in the present paper (see e.g. inset of Fig. 5). This experimental 
study also shows that the particle transit time is Gaussianly distributed in time (see Fig. 9 in ||]) in accordance with 
the solution @ of Eq.(jTf) (see Fig.l). 

Front-type dynamics is also encountered in shock propagation in homogeneous or inhomogeneous media fill or 
packet transport in the Internet [ fl2| . As the propagation-dispersion equation is for the first-passage time distribution, 
it should also be suited for the description of transport driven by an input current in a disordered random medium 
p3[ . In the area of traffic flow, there are typical situations where cars moving on a highway from location A to location 
B, are subject to time delays along the way, and - with the assumption that all cars arrive at destination - one wants 
to evaluate the time of arrival [Q . Financial series as in the time evolution of stock values are another example [|l5| : 
over long periods of time (typically years) one observes a definite trend of increase of, for instance, the value of the 
dollar. So any preset reachable value will necessarily be attained, the questions being: when? and what is the time 
distribution around the average time for the preset value? While the classical question is: after such or such period 
of time, which value can one expect?, there might be instances where the reciprocal question should be considered. 
Because of the generality of the propagation-dispersion equation, it should be expected that, either in its simple form 
( |l8| ) or in its generalized form (p2|), the equation will be applicable to a large class of first-passage type problems in 
physics and related domains. 



We thank Alberto Suarez and Didier Sornette for useful comments. JPB acknowledges support by the Fonds 
National de la Recherche Scientifique (FNRS, Belgium). 



The t^s in the first visit equation, Eq. (JT|) , are known explicitly for specific microscopic dynamics of a particle on a 
lattice 0. Two parameters are used to characterize the trajectory of the particle in the propagation channel. First, 
as the propagation line is not necessarily along one of the lattice axes, we define m (> 1) as the minimum number 
of time steps required to travel from site I to the neighboring site I + 1, i.e. to perform a displacement Sr along 
the propagation line. The second parameter, a, is the length of an 'elementary loop', i.e. the minimum number 
of displacements necessary to return to a site. Typical values of m and a are shown in Fig. 4 for various channel 
geometries. The expression for the time delays then reads fj 



For instance, in the one-dimensional lattice (where n = 1, a = 2, and m = 1; see Fig. 4a) |5J] j ' = corresponds 
to straightforward motion from site I to site I + 1 in one time step, and J = 1 corresponds to one step backwards 
followed by two forward steps. In the square lattice, a — 4, and m — 8 (see Fig. 4c) for Langton's ant dynamics j|, [|. 
The value m — 8 follows from the fact that when the particle arrives at site 1 for the first time, the shortest possible 
path to the next first visited site on the edge of the propagation channel goes to site 2 as shown in Fig. 4c. Indeed, for 
Langton's ant dynamics, all sites are initially in the same state (here scattering of the particle to the left). So when 
the particle visits site 1 for the first time, the site located immediately North of site 1 has not yet been visited and 
is therefore in the left-scattering state. It is then easy to figure out that for any path different from the path shown 
with heavy solid line in Fig. 4c, the particle will make a longer excursion to arrive at site 2. 

With the specification of the time delays in (|55| ) , we can give an explicit formulation of the quantities which appear 
in the continuous limit equation. The average displacement time and the variance read respectively 



Acknowledgments 



Appendix: Specific microscopic dynamics 



(55) 



1 1 




(56) 



3=0 



and 
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a = 2 m = 1 



a) 



b) 




1 2 3 4 5 6 



c) 




FIG. 4: Propagation channels in one- and two-dimensional lattices. Heavy lines with arrows show particle displacements in 
elementary loop (a — 2, 3, 4) and minimal trajectory between two consecutive sites on the propagation line (m = 1, 2, 8). 
Propagation in the one-dimensional lattice (a) and in the triangular lattice (b) is discussed in j^]; (c) refers to Langton's ant 
(see || and |B|). 
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where from it follows that 



and 



J=0 3=0 J 

= ({j 2 )-{j) 2 )a 2 m 2 (5t) 2 , (57) 



- = (i + a (j»m£, (58) 
c or 



7=((i 2 )-(j) 2 )(«m) 2 M. (59 ) 

or 

The value of c and of 7 depends on the spin orientation probability q, the probability that a site be in the | state at 



the initial time. The derivation of the propagation-dispersion equation proceeds exactly along the lines of Section IV 
and yields Eq. (|l^) with c and 7 given by ( |58| ) and (^9|) respectively. Figure 5 shows the agreement of the analytical 
results with the simulation data for the 1-D and 2-D lattices. For instance, for the 1-D spin-lattice whose dynamics 
is described in the introductory section, m = 1, a = 2 (see Fig. 4a), j — 0,1, and (j) — J2j3Pj = 1 ~ <Z an d 
(j 2 ) — (j) 2 = J2jJ 2 Pj ~ (J2j JPj) 2 = ^(l — so c = (3 - 2g) _1 and 7 = 4q(l — q). The corresponding numerical 
values are given in the figure caption for both the 1-D lattice (Fig. 5a) and the 2-D triangular lattice (Fig. 5b). 

In a recent paper (l6| , Bunimovich and Khlabystova, referring to a preliminary version of the present work [ fl7| , 
obtain the same propagation-dispersion equation for the case of a particle moving in a rigid environment |Ts| : the 
rigidity factor f is defined as the number of visits of the particle to a site necessary to flip its state (T^=H)- The 
time delay probabilities are then space-dependent. For odd rigidity (NOS model in H), there is no qualitative 
difference for c and 7 with the simple 1-D case (f = 1) discussed above (compare the expressions given in the previous 
paragraph with those in Section 4.1 of Jig] ). However when the rigidity factor takes even values, the propagation 
speed and the dispersion coefficient become space-dependent (see Section 4.2 in Jl(|). So the microscopic dynamics 
of the (NOS) model with even rigidity offers an example where the continuous limit is given by the generalized 
propagation-dispersion equation (p2|). 

The explicit expression ( pS] ) for Tj based on the automaton dynamics yields analytical expressions for c and 7 
in terms of the spin-lattice characteristics and there from in terms of the probability q. Concomitantly there is an 
explicit reference to a feed-back mechanism where the dynamics modifies locally the substrate which in turn modifies 
the dynamics, and there are systems where this mechanism should be important. However we emphasize that the 
specification ( |55| ) of the delay time Tj is not indispensable. The propagation-dispersion equation (jl8) is general as it 
follows from the continuous limit of the first visit equation (jj]) without recourse to ([55]) as shown in Section IV. It 
suffices that there exists a distribution of time delays with X^=o Pj ~ 1 anc ^ with finite moments (r a ), to obtain ( p^ ) 
from (Q) . This establishes the validity of the propagation-dispersion equation for a class of systems whose dynamics 
is subject to time-delays independently of the underlying microscopic mechanism responsible for the delays. 
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